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In this paper we introduce a new method to design interparticle interactions to target arbitrary 
crystal structures via the process of self-assembly. We show that it is possible to exploit the curvature 
of the crystal nucleation free-energy barrier to sample and select optimal interparticle interactions 
for self-assembly into a desired structure. We apply this method to find interactions to target 
two simple crystal structures: a crystal with simple cubic symmetry and a two-dimensional plane 
with square symmetry embedded in a three-dimensional space. Finally, we discuss the potential 
and limits of our method and propose a general model by which a functionally infinite number of 
different interaction geometries may be constructed and to which our reverse self-assembly method 
could in principle be applied. 



INTRODUCTION 

Not only is understanding, controlling and predicting 
the phenomenological behavior of particle self-assembly 
one of the great mathematical challenges for the 21st 
century, but its applications in materials science and en- 
gineering hold promise for the development of materials 
with novel electronic, mechanical, and optical properties. 
Although most of the work performed in this field is his- 
torically rooted in the self-assembly of small molecules, 
the last decade has witnessed extraordinary advances in 
particle synthesis at the meso-scale [TH5] , making possible 
the production of building blocks with complex chemical 
and geometrical properties with an unprecedented degree 
of precision. Unfortunately, a coherent theoretical frame- 
work around the problem of self-assembly is still missing, 
and numerical simulations have taken the lead in explor- 
ing the wealth of new phenomenological behavior arising 
from the collective behavior of non-isotropic components. 

Most numerical studies on self-assembly of nanopar- 
ticles performed so far have adopted the patchy sphere 
model [6 . In this model, the isotropicity of a particle is 
broken by placing on its surface regions (patches) with 
different physical properties; for example, hydrophobic 
chemical groups or single- stranded DNA chains. Theo- 
retically, these regions are incorporated into the inter- 
particle potential by a simple angular dependence which 
favors or disfavors the alignment of such patches. Al- 
though self-assembly of several simple structures has 
been achieved with the patchy models (references [7]- 
[9] are just a few examples of the large body of work 
published on the subject; for a recent review see refer- 
ences [l0l[TT]), a general modeling approach to the prob- 
lem is missing. Shape and position of the interaction sites 
is either guessed using physical arguments, or inspired by 
known molecules or protein structures aggregating into 
a similar target crystal. There are two notable excep- 
tions: the inverse optimization technique proposed by 
Torquato [12], which is specific to nondirectional inter- 
actions, and so-called "bottom-up building block assem- 



bly," devised by Jankowski and Glotzer [13], which re- 
quires the construction of the most relevant terms of the 
partition function of the system, starting from individual 
particles. The development of an efficient numerical pro- 
cedure to design interactions between nanoparticles that 
targets specific crystal structures via the process of self- 
assembly would therefore be a result of great importance. 

Although the generic features of particle aggregation 
can be described, at least phenomenologically, in terms 
of simple thermodynamic arguments [T4HT8] , the details 
of the self-assembly process are far from being under- 
stood, even in the simple case of aggregation of isotropic 
particles into macroscopic three-dimensional crystals. In 
fact, a full theoretical description of this problem must 
incorporate critical kinetic effects which are not captured 
by classical thermodynamics, and which have dramatic 
macroscopic consequences [9j [19] [20] . 

It is now understood that for self-assembly to take 
place, a very delicate balance between entropic and ener- 
getic contributions, coupled to a precise geometric char- 
acter of the components, must be satisfied. In general, 
self-assembly of nanocomponents is not to be expected 
unless a careful design of the building blocks has been 
performed beforehand. [9] [19] 

Figure [T] illustrates the problem for a simple model: 
spherical particles with attractive patches oriented to 
form a two dimensional honeycomb lattice. When the 
angular size of each patch is too large, the interaction 
is not specific enough to select the desired crystal, result- 
ing in amorphous structures originating from the compe- 
tition of multiple fitting geometries. If Q is too small (too 
specific) , the probability that two particles in close prox- 
imity are properly aligned to interact becomes negligible, 
and the system is found in a gas phase unless a very large 
interparticle energy E is provided, which in turn drives 
the system into a gel phase. Analogous arguments can 
be made for the overall strength of the interaction that, 
for a reasonable patch size, should be neither too strong 
nor too weak. The net result is that self-assembly is a 
very elusive process that requires a careful design and 
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FIG. 1: Illustration of a generic self-assembly diagram for 
patchy spherical particles expected to aggregate into a hon- 
eycomb lattice. E is the interparticle attraction strength and 
Q is the angular size of the patches. 



fine tuning of the interparticle interactions, and typically 
the target region in the interaction space in quite narrow. 

Our interest is in the problem of "reverse self- 
assembly," in analogy to the problem of "reverse pro- 
tein folding" in which protein sequences are designed to 
yield a desired ground-state folded structure [2lJ[22]. The 
problem can be formulated as follows: given an arbitrary 
final structure, is there an efficient method by which in- 
terparticle interactions may be designed so that the parti- 
cles spontaneously self-assemble into the given structure? 

In this paper we use simple physical arguments to 
develop a numerical procedure capable of sampling the 
space of interactions, in terms of both patch geometry 
and binding energy, to generate nanoparticle interactions 
leading to self-assembly of desired crystal structures. 
Without loss of generality, we limit our discussion to 
spherical particles interacting via anisotropic short-range 
attractive pair potentials mimicking the hydrophobic in- 
teractions driving self-assembly of Janus particles [9l [23] . 

Classical nucleation theory provides a simple frame- 
work within which to think about crystal formation. A 
free energy gain (/i c — Pf)n is associated with the for- 
mation of a nucleus of n particles of the crystal phase 
at chemical potential \i c out of a fluid phase at a larger 
chemical potential pf. A free energy cost jA is associ- 
ated to the formation of an interface of area A (ex n 2 / 3 ) 
and surface tension 7 between the two phases. Minimiza- 
tion the total free energy AG with respect to n leads, for 

a spherical nucleus, to a critical size n c = (^j^j^J (p 
is the equilibrium crystal density). Crystal nuclei larger 
than n c will grow until the phase transformation is com- 
pleted, all others will shrink and vanish. We argue that 



a successful strategy for crystal design should take into 
account the physical properties of the parent fluid phase, 
and our working hypothesis is that the free energy of 
crystal nucleation can be exploited to design interactions 
to target arbitrary crystal structures. The main idea is 
to force a crystal nucleus of desired symmetry to be in 
contact with its own fluid and use a numerical procedure 
to select for those interactions between the particles that 
minimize the free energy cost required to hold that nu- 
cleus in place. Our scheme consists of two parts: 1) we 
determine the optimal shape of the hydrophobic regions 
(ft) satisfying the condition stated above, and 2) given Q, 
we find the interaction strength (E) for which the system 
is likely to nucleate into the target (defect-free) crystal 
phase. 

SAMPLING THE INTERACTION SHAPE Q 

Consider a system of N identical particles with a given 
interparticle potential U (O^, r) set in a volume V. Define 
an order parameter q capable of detecting the symmetry 
of the desired crystal phase. Grow from the fluid and 
equilibrate a crystalline nucleus of size no using a stan- 
dard bias Monte Carlo method targeting the size of the 
largest crystalline cluster in the system, n, via a poten- 
tial V B (n) = f (n - n ) 2 [24 . Set the binding energy 
among the particles to a sufficiently small value to en- 
sure that the nucleus melts once the bias is removed, and 
compute from a full simulation in the presence of the bias 
the average crystal size n(Qi). 

Now define a design potential Vb[n(f2i)] = — cm(^), 
where a is a numerical constant. At this point the idea 
is to sample over the space of interactions using Vd as 
a driving force. Specifically, we generate an alternative 
(trial) shape for the interaction between any two particles 
in the system Qj = Qi+AQ and repeat the previous steps 
to obtain a new estimate for Vo[n(^j)]. Qj is accepted or 
rejected based on a standard Metropolis criterion, thus 
ensuring that the Q will be driven towards values that 
maximize the size n, i.e. minimize the load requested of 
the bias to hold the crystal in place. 

SAMPLING THE INTERACTION STRENGTH E 

Unfortunately, our method does not allow easy mea- 
surement of the height of the nucleation barrier given 
an interaction strength E. This is mostly because the 
surface tension between the crystal and the fluid phase 
is unknown. Nevertheless, we have direct access to the 
slope of the free energy barrier. Therefore, although one 
cannot design the system to comply with a specific nu- 
cleation rate z/, by modulating E one can design the size 
of the critical nucleus n c . A critical nucleus that is too 
small will result in the almost instantaneous nucleation 
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of several crystallites that will form defects and grain 
boundaries as they meet while growing. The opposite 
scenario will lead to absence of crystallization within the 
experimental time frame. For the systems we have ex- 
amined, we find that n c ~ 15 — 30 results in nucleation 
events that are quick, yet sufficiently rare to prevent for- 
mation of multiple crystals. The choice of n c may require 
a few iterations depending on the details and the size of 
the system. 

The strategy behind the design of the critical nucleus 
size is analogous to that described in the previous case, 
except that the design potential in this case has a har- 
monic functional form defined as Vb[n(^)] = a(n(Ei) — 
n c ) 2 , and we sample over the interaction strength E. 
Minimizing Vd implies that the system will be driven 
towards that value of E (E c ) for which the nucleus has 
the same probability of growing or shrinking. This con- 
dition guarantees that the system is at the top of the 
nucleation free energy barrier, and that n c has indeed 
become the critical nucleus by definition. 

Note that in principle, the interaction geometry could 
have an arbitrary number of parameters that could all be 
optimized simultaneously; however, it is crucial that the 
optimization of the geometry precedes the optimization 
of the strength of the potential. In fact, it is mandatory 
for the nucleus to be precritical in order for the geometry 
optimization scheme to be effective. 



NUMERICAL TESTS 

As a proof of concept for our method we consider the 
design of two distinct crystal structures for which we can 
guess the solution in the interaction space and know how 
to define an order parameter q: a simple cubic crystal 
(SC) and a two dimensional sheet with square symmetry 
embedded in a three dimensional environment (2SQ). 

For both systems we adopt the Kern-Frenkel model [6 . 
Particles are described as hard spheres of diameter a in- 
teracting with a short-range attractive interaction that is 
turned on whenever hydrophobic regions (the patches) on 
different particles face each other. For each pair of par- 
ticles i and j with patch indices a and /3, the interaction 
is defined as 



(1) 



where ^sw(jr/) is a standard attractive isotropic square 
well potential of depth e and range 1.15cr, and f ap (Qij) 
depends on the particles' mutual orientations and is de- 
fined as 



i if 



iij • e a > cos 
and Yji • e# > cos# (2) 



Here 6 is the angular size of the hydrophobic regions 
(selected to be all identical in size and circular in shape), 
ry is the unit vector along the direction of the interpar- 
ticle separation, and e a is the unit vector connecting the 
center of a particle to the center of the patch a on its 
surface. 

In these simple systems and e are the design param- 
eter we need to tune for self-assembly to take place. All 
of our simulations are performed in the NVT ensemble 
using a minimum of 256 particles in a box with periodic 
boundary conditions. A good order parameter to detect 
SC crystals is the standard local bond order based on 
spherical harmonics, [24-26 . Given a particle z, we 
compute 



1 N b {i) 



(3) 



where j runs over the Nb(i) neighbors of particle z, from 
which a rotationally invariant order parameter correlat- 
ing the orientation of neighboring particles i and j can 
be defined as 



(4) 



Once averaged over all neighbors j, the resulting quan- 
tity (?4 is compared to a cutoff, g cut , to decide whether a 
particle can be tagged as crystalline or not. 

For the 2SQ case we used q^ with the added constraint 
that a particle must have interactions with no more than 
four neighbors in order to be considered "crystalline." 
The location of the patches automatically prevents the 
formation of SC crystals in this case. The insets in Figure 
2 sketch the patch positions over the particles. 

Figures [2 (a) | and |2(b)| illustrate how the force F(0) = 
— n{n — no) required to hold a nucleus of no particles 
immersed in its fluid phase depends on the size of the 
circular regions for the SC and the 2SQ crystals re- 
spectively, and specifically Fig. |2(a)| also shows that the 
optimal value is fairly independent of the particular size 
of the nucleus no held in contact with the fluid. 

The corresponding simulations were performed at den- 
sities of psc = 0.1 and p2SQ = 0.01, binding strength 
£sc = 3.5/cbT and E2SQ — 5.75&b^ 5 and a harmonic 
bias potential of spring constant ksc — 0.2/cbT and 
fesQ = 0.4/cbT. Clearly, F(&) is a sufficiently sensitive 
parameter to discriminate among the different angular 
sizes, and presents in both cases a distinct optimal value; 

also shows that 



J sc 



22 c 



and 6* 2SQ 



20°. Figure 2(a) 



otherwise 



the optimal value is fairly independent of the particular 
average size no of the nucleus held in contact with the 
fluid. Figure [3] shows how the location of 0* and £* can be 
obtained automatically by using the Monte Carlo scheme 
in the space of interactions. 
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(a) 



(b) 



FIG. 2: Force vs. for (a) SC and |(b)| 2SQ cry stals . The insets show snapshots of the target crystals, and sketches of the 
locations of the patches in our particle model. In (a) the different lines represent data obtained by imposing different nucleus 
sizes no as indicated in the legend. 
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FIG. 3: Monte Carlo trajectories in the space of interactions for the design of the simple cubic crystal. In (a) the shape of the 
patches defined by the solid angle is allowed to fluctuate while keeping the strength of the interaction e constant. In (b) e 
fluctuates while keeping constant and at the optimal value found in (a). 



It should be stressed that the minimization of 
Vd[^(^)] can b e achieved using any minimization algo- 
rithm; nevertheless, we find that the Monte Carlo scheme 
allows us to use shorter simulations, for each trial 0$, 
than what would be necessary for other direct minimiza- 
tion schemes. The reason is related to the precision of 
the estimate of n for relatively short trajectories that 
could be over- or underestimated. This could lead to 
fictitious local minima, which could trap a direct mini- 
mization scheme, but are easily overcome with a standard 
Monte Carlo method. 

In order to check our method, e — phase dia- 



grams were constructed using the traditional, "forward" 
method of trial-and-error, running Monte Carlo simu- 
lations for 10 7 steps and determining whether crystal- 
lization occured. As shown in Figure [4j the parameters 
detected by our methods are within the crystallization 
regions for the two target crystals. Unsurprisingly, the 
result falls roughly to the high-0, low-e edge of the crys- 
tallization region; recall that was selected using a value 
of e too low for crystallization, which would be expected 
to result in a larger 6 (note the roughly inverse e — 6 re- 
lationship in Figure [4]) ; e was then selected using the 
found in the first step. 
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(a) 



(b) 



FIG. 4: Phase diagrams for (a) SC and |(b)| 2SQ crystals. Lines show the border around the phase region in which particle form 
the desired crystal; the points marked by an 'X' are the (0, /3s) combinations found to be optimal by our method. 



LIMITATIONS 

It is important to discuss the limitations of the method. 
First of all, in its actual formulation, our method only 
works for systems that will self-assemble into an infinitely 
large aggregate via the process of nucleation. It is not 
obvious how to generalize it to include self-assembly into 
finite size aggregates such as for instance viral capsids. 

The second limitation is that although the method pro- 
vides a solution to the reverse self-assembly problem, 
there is no guarantee that the solution is the optimal 
one. This is because our method forces the nucleation 
process to follow the classical route, i.e. the forming nu- 
cleus has the same structure as the target crystal; how- 
ever, there are several examples [H [20j I27H29] where the 
nucleation barrier may be lowered by following a more 
complex dynamical pathway that may include metastable 
states having different symmetry than that of the target 
crystal. For instance, it is possible to imagine that the 
formation of the SC crystal could benefit from an addi- 
tional weak, non-specific, isotropic interaction, on top of 
that provided by the patches, that may initially lead the 
system into a high-density metastable fluid phase from 
which nucleation into the final structure may proceed at 
a faster rate than that predicted otherwise. 

Finally, it is crucial to develop a good order parameter 
q to describe the desired crystal structure. Figure [5] il- 
lustrates how an inefficient order parameter may lead to 
fictitious minimization in the space of interactions while 
designing the 2SQ crystal. The different lines in the F 
vs. diagram represent different values of q cut (denning 
how restrictive the order parameter is) from 0.8 to 0.99. 
We find that a cutoff in the order parameter of at least 
0.97 is required to adequately distinguish between the 
square and hexagonal symmetries for large values of 0. 




FIG. 5: F vs. diagrams. The dependence of the method 
on the order parameter used. Different curves correspond to 
different values of the cutoff q cu t- 



The curves related to the less restrictive order parameters 
would in fact misleadingly indicate a flatter bottom of the 
curve, while in reality we find that any angle larger than 
~ 25° will lead to nucleation into a two-dimensional crys- 
tal with hexagonal symmetry. Adding an energy penalty 
to prevent arrangements compatible with the competing 
six- fold symmetry (apart from imposing a limit to the 
number of neighbors) may also be considered as a means 
of improving the design procedure. 



BEYOND THE KERN-FRENKEL MODEL 

Here we propose a more general model to describe in- 
terparticle interactions that we name the Adaptive Pixel 
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Model The idea stems from the need to devise a way of 
sampling over different geometric patterns (beyond cir- 
cular) in search of those which can efficiently hold the 
single components into a desired target structure. The 
first step is the discretization of the surface of the parti- 
cle. 

For spherical building blocks, we cover the surface 
of each particle with a large number, N p , of regularly- 
spaced interaction sites (pixels) as illustrated in Fig. [6] 
A good arrangement of the pixels can be obtained by 
using the spherical triangulation provided by an (n, m) 
delta-icosahedron [30 , and N p is selected depending on 
the complexity of the target structure. Euler's theo- 
rem [30 imposes the following constraint on N p : N p — 



10(n 2 



m 2 ) + 2 [30 , where n and m indicate that 



one has to move n pixels along the row of neighboring 
bonds on the sphere, and then after a turn of 120°, move 
for m extra pixels. 

In the simplest version of the model, to each pixel k on 
a particle z, is assigned a variable which has a binary 
character, G {1,0} depending on whether that inter- 
action site is switched on or off. Whenever two particles 
i and j are within a given distance of each other, the axis 
between them, r^-, is calculated. If the nearest digit to 
the point where nj crosses the surface of each particle 
is on, then the particles feel an overall short-range at- 
tractive interaction. Pixels on the same particle do not 
interact with each other. The interaction pair potential 
between any two particles, i and j, of diameter a, set at 
a distance nj from each other, then takes the form 



V{n 3 







if < r 
otherwise 



(5) 



where and Sji are the binary variables corresponding 
to the digits intersected by nj on particles i and j, re- 
spectively, as described above. Excluded volume between 
the particles is enforced via a standard hard-sphere po- 
tential, Vrs- 

The main advantage of this setup is that once parti- 
cles are held into place at given positions, the geometry 
of the interacting regions emerges as a result of a simple 
Monte Carlo simulation on which samples different 
states according to Eq. [5] Crucial to the efficiency of the 
model, is the independence of the interaction strength 
of the total area of the attracting region. This condi- 
tion, also assumed in the Kern-Frenkel model, is appro- 
priate when considering interactions that have a range 
of action that is small compared to the colloidal diame- 
ter ((r*o — a) i$ 0.15a), and allows us to circumvent the 
overwhelming cost related to the computation of the TV 2 
distances between the pixels. Furthermore, as the rel- 
ative distances of on-particle pixels are frozen and only 
active pixels need to be tracked, it is possible to perform 
the search of the nearest pixel to any point on the sphere 
very efficiently. This is achieved by creating a cell list 




FIG. 6: Illustration of the structure of the Adaptive Pixel 
Model, on pixels are depicted in red while off pixels are in 
gray. The magnification in the top image shows the Voronoi 
tessellation around the pixels (computed as described in the 
text). The effective geometry of the active sites in this repre- 
sentation is a hexagon. 



over the spherical particle surface in and (f) (the spheri- 
cal coordinates), and by associating to each cell the iden- 
tity of the nearest pixel. This is equivalent to generating 
a discrete Voronoi tessellation [30] of the spherical sur- 
face based on the pixel locations (see sketch in Figure [6]), 
which needs to be performed only once at the beginning 
of the simulation. Any shape for the interaction regions 
can be achieved by simply switching on or off pixels or 
groups of pixels. 



CONCLUSIONS 

In this paper we have proposed a simple two-step 
method for the problem of reverse self-assembly. The 
idea is to exploit the curvature of the nucleation free- 
energy barrier to sample and select optimal inter-particle 
interactions for self-assembly into a target structure. We 
presented numerical simulations to test the efficacy of 
our method, and discussed in detail its limitations and 
its potential. These simulations show that our method 
reduces the time to solve the problem of determining op- 
timal interaction parameter from on the order of weeks 
(for trial-and-error approaches) to hours. 

Finally, we proposed a new model, the Adaptive Pixel 
Model, by which almost any interaction geometry can 
be realized in a simple and efficient way. It should be 
stressed that our method is not limited to spherical par- 
ticles but can be applied to any particle shape. In princi- 
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pie, particle shape could be introduced as a new param- 
eter in the interaction space and be sampled over using 
the scheme proposed in this paper. 

Although our method does not capture the dynamical 
subtleties of the crystal formation process, it does provide 
a very efficient way of screening over a large number of 
given interaction geometries that can be mapped onto the 
pixels. Efficient ways of sampling the interaction space 
could be obtained using genetic algorithms that can be 
used to evolve optimal interaction patterns given a set of 
initial shapes. Work in this direction is currently under 
investigation and will be published elsewhere. 
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